sor_solver_linearsys <- function(A, b) {
  x <- matrix(1, nrow = ncol(A))
  x_new <- x
  for (j in 1:72) {
    for (i in 1:ncol(A)) {
      x_new[i] <- x[i] + (b[i] - A[i,]%*%x_new)/A[i, i]
    }
    x <- x_new
    print(x_new)
  }
  return(x_new)
}

A <- matrix(c(4, -2, -1, -2, 4, -2, -1, -2, 3), 3, 3)
b <- matrix(c(0, -2, 3), 3, 1)

solve(A, b)


x <- sor_solver_linearsys(A, b)

sor_solver(A, b, matrix(-Inf, 3), matrix(Inf, 3), matrix(1, 3), max.steps = 2)


conjugate_gradient_method(A, b, matrix(1, 3), 4000)
